UNCLASSIFIED 


AD  NUMBER 

AD462335 

NEW  LIMITATION  CHANGE 
TO 

Approved  for  public  release,  distribution 
unlimited 


FROM 

Distribution  authorized  to  U.S.  Gov't, 
agencies  and  their  contractors; 
Administrative/Operational  Use;  FEB  1965. 
Other  requests  shall  be  referred  to  Bureau 
of  Naval  Weapons,  Washington,  DC. 


AUTHORITY 

NAVWEPS  ltr,  2  Feb  1961 


THIS  PAGE  IS  UNCLASSIFIED 


m 

CO 

00 

i^TG-648 
^EBRUARY  1965 
Copy  No.  a 


CD 


O  r-s 
LU  S~l 
CD  52 

S  C3 

«=C  CSC 

S  CO 

CJ>  *=C 


•*»? 


cbt*» 


■Si?Su 

ddc  *  **  - 


Technical  Memorandum 

A  FORTRAN  COMPUTER  PROGRAM 
FOR  THE  SOLUTION 
OF  MULT^HMENSIONAL. 


TRANSECT,  ABLMiGN  PROBLEMS 


by  R.  W.  ALLFN  and  R.  P.  SUESS 


D&C 


it 


§ 


THE  JOHNS  HOPKINS  UNIVERSITY  •  APPLIEw  PHYSICS  LABORATORY 


Best 

Available 

Copy 


TG— 648 

FEBRUARY  1965 


Technical  Memorandum 

A  FORTRAN  COMPUTER  PROGRAM 
FOR  THE  SOLUTION 
OF  MULTI  DIMENSIONAL, 
TRANSIENT,  ABLATiON  PROBLEMS 

by  R.  W.  ALLEN  and  R.  P.  SUESS 


THE  JOHNS  HOPKINS  UNIVERSITY  •  APPLIED  PHYSICS  LABORATORY 
8621  Georgia  Avenue,  Silver  Spring  Maryland  20910 

unde  Control  NOw  64 -0604 -c,  Bvrew  of  Weepom,  Depart--  of  the  Nvvy 


i- 


ABSTRACT 

A  flexible  FORTRAN  computer  program  to  d.-ccimlne  the  temperature  history 
and  ablation  history  of  aerodynamically  heated  flight  bodies  has  been  devised. 

The  effect  on  heat  conduction  and  aerodynamic  heating  histories  of  the  removal 
of  ablating  wall  materials  is  automatically  accounted  for  by  the  computer  pro¬ 
gram.  The  program  is  flexible  in  that  it  can  accommodate  flight  body  heat  con-- 
ductiou  inputs  in  the  form  of  any  reasonable  combination  of  geometry  and 
construction  materials.  The  program  will  accept  only  one  ablative  material  at 
any  single  position  on  the  flight  oody  surface  but  different  ablative  materials 
can  be  specified  for  different  locations.  One  section  of  the  program  receives 
flight  trajectory,  radiation-property ,  and  local  aerodynamic  flow  inputs  end 
provides  for  the  computation  of  local  aerodynamic  heating  ar.d  radiation  relief. 
The  other  section  of  the  program  governs  the  computation  of  temperature  history 
throughout  the  flight  body  and  the  thickness  history  of  ablating  layers.  There 
is  no  provision  for  readjusting  vehicle  aerodynamics  in  accordance  with  body- 
shape  changes.  In  setting  up  the  program  it  was  assumed  that  the  process  of 
decomposition  of  the  ablative  material  was  concentrated  at  the  surface  and  that 
the  chemistry  of  the  decomposition  could  be  accounted  for  by  empirically-based 
effective  heats  of  ablation  and  ablation  surface  temperatures  instead  of  chemical 
reaction  equations  and  associated  chemical  reaction  rate  constants.  Since  the 
report  is  to  serve  as  a  user's  manual.,  a  detailed  description  of  the  computer 
program  is  provided.  The  description  covers  the  engineering  »,•£''  o»,«h  adopted  in 
relation  to  the  general  problem  of  determining  temperature  e  id  ablation  historiea 
in  a  flight  body  by  numerical  methods.  This  is  foixowed  by  -information  on  the 
program  structure,  flow  of  program  information,  and  P0RTRAN  nomenclature .  A 
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samp le  problem  is  posed  and  step-by-step  preparation  of  the  F0RTRAN  codi) 
is  explained.  Finally,  an  actual  print-out  from  a  computer  run  for  the  : 
problem  is  dis" '  r.  ed  and  discussed. 
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A  F0RTRAN  COMPUTER  PROGRAM  FOR  THE  SOLUTION  OF 
MULTI-DIMENSIONAL,  TRANSIENT,  ABLATION  PROBLEMS 

INTRODUCTION 


One  of  the  many  considerations  in  the  design  of  a  hypersonic  defense 
missile  is  the  thermal  protection  system.  Ablative  systems  are  attractive  in 
this  regard  because  of  their  automatic  response  to  aerodynamic  heat  flux  and 
consequent  simplicity  compared  to  transpiration  or  ducted  cooling  systems.  As 
a  result  of  this  situation,  there  is  a  continuing  demand  for  efficient,  flexible, 
and  effective  ablation  design  calculation  procedures.  This  demand  has  been  met, 
in  part,  by  the  appearance  of  a  number  of  ablation  computer  programs.  However, 
such  programs  have  tended  to  be  oriented  toward  the  re-entry  problem  and  it  h.is 
seemed  that  more  attention  could  have  been  given  to  constructing  programs  flexi¬ 
ble  enough  to  accept  the  variety  of  geometries,  materials,  and  trajectory  inputs 
associated  with  the  hypersonic  missile  problem.  This  report  presents  an  ablation 
computer  program  for  subliming  ablators  that  has  been  constructed  with  this 
flexibility  in  mind.  The  report  was  prepared  both  as  a  descriptive  exposition 
of  a  particular  engineering  method  and  as  a  user's  manual.  Thus  it  lays  out  for 
the  potential  user,  not  only  the  general  engineering  techniques,  but  also  the 
flow  of  program  information  and  the  step-by-step  instructions  to  the  person 
preparing  the  F#RTRAN  coding  sheets.  Before  proceeding  further,  a  quick  review 
of  ablation  processes  is  in  order. 

As  is  well  known,  ablation  is  characterized  by  the  disintegration  of  s 
solid  material  under  the  combined  mechanical  and  thermal  action  of  a  boundary 
layer  of  hot  gas.  In  high-speed  flight,  ablation  generally  occurs  as  a  ~.om- 
bination  of  physical  and  chemical  processes.  The  principal  physical  modes  of 


ablation  are: 
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(a)  Melting  of  a  solid  surface 

(b)  Vaporization  of  a  liquid 

(cl  Sublimation  of  a  -=■  >  1  i d  surface 
(dl  Fra  curing  or  shearing  of  solid  surface  layers 
Chemical  processes  occurring  during  ablution  fsll  into  the  following  major 
categories : 

(a)  Oxidation  of  solid  surface 

(b)  Oxidation  of  released  gas 

(c)  Molecular  break-down  of  solid  under  action  of  heat  (pyrolysis) 

In  typical  ablation  processes  various  combinations  of  the  foregoing  reactions 
usually  occur.  Some  examples  are  the  following: 

Quartz  -  The  solid  surface  molts  and  is  swept  along  by  the  air  boundary 
layer  while  the  liquid  vaporizes  into  the  boundary  layer. 

Teflon  -  The  supermolecule  (polymer)  breaks  down  in  depth,  releasing 
monomer  gas  and  the  gas  is  oxidized  at  the  surface  by  oxygen  from  the 
boundary  layer. 

Graph: ce  -  Oxygen  diffuses  inward  through  the  bound  try  layer  to  support 
combustion  of  carbon  at  the  graphite  surface. 

Phenolic  Resin  -  'Tie  supermolecule  breaks  down  in  a  high -temperature 
surface  layer.  Gases  are  released  leaving  a  porous  carbon  (char)  layer. 
The  released  gases  react  with  oxygen  which  is  diffusing  inward  through 
the  air  boundary  layer. 

In  designing  the  computer  program,  attention  was  focused  on  the  case  of  the 
pyrolyzing  ablator.  It  was  assumed  that  volume-loss  calculations  would  be 
initially  restricted  to  the  case  in  which  the  exposed  surface  of  the  ablator 
undergoes  a  pyrolyzing  process  equivalent  to  a  sublimation  phase-change.  It  was 
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anticipated  that  more  complex  ablation  processes  such  as  charring  could  be 
accommodated  at  a  later  date.  Meanwhile,  however,  if  surface  combustion  takes 
place  cr  if  chemical  decomposition  takes  place  in  depth,  the  sicuaticn  is  accom¬ 
modated  by  assuming  the  chemical  decomposition  processes  to  be  concentrated  at 
the  surface  in  the  form  of  a  sublimation  phase-change  having  a  known  effective 
heat  of  ablation  and  a  known  ablation  surface  temperature.  As  a  result,  the 
computer  program  is  realistically  based  on  the  availability,  in  the  literature, 
of  data  on  the  effective  heat  of  ablation  and  ablation  surface  temperature. 

The  numerical  technique  is  based  on  the  familiar  lumped -parameter  approach 
in  forward-difference  form  and  is  programmed  specifically  for  the  IBM  7094  digital 
computer.  The  program  consists  of  a  number  of  FORTRAN  and  FAP  subroutines  which 
are  called  into  action  upon  command  of  a  FORTRAN  control  program.  The  control 
program  is  prepared  by  the  user  in  accordance  with  the  problem  at  hand.  In  a 
broad  sense,  the  program  combines  boundary -layer  local  heat  transfer,  ablation, 
and  radiation  at  the  surface  with  multi-dimensional  heat  conduction  beneath  the 
surface  of  a  body  in  non-steady  high  speed  flight.  The  program  accounts  for  the 
effect  of  in-flight  body  contour  changes  on  non-steady  heat  conduction  but  does 
not  account  for  the  effect  of  in-flight  body  contour  changes  on  vehicle  aero¬ 
dynamics.  Hence,  addition  of  this  feature  is  planned  for  sometime  in  the  future. 
Input  data  consist  of  flight  trajectory  and  flight-body  details  and,  in  particular, 
the  ablation  temperature  and  a  modified  form  of  the  effective  heat  of  ablation. 

The  report  is  divided  into  four  differed  sections  ranging  in  content  from 
engineering  considerations  to  a  solution  of  an  actual  ablation  problem.  The 
organizational  plan  is  as  follows : 

Section  1.  Problem  Formulation  and  Governing  Equations.  This  section 
covers  the  engineering  methods  employed. 
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Section  2. 


Section  3. 


Section  4. 


Structure  of  Program.  This  section  covers  the  inputs,  sub¬ 
programs,  ana  information  flow. 

The  F0RTRAN  Control  Program.  This  section  describes  the  content 
of  the  principal  F0RTRAN  statements  and  gives  the  program 
nomenclature . 

Sample  Problem.  This  section  describes  the  detailed  preparation 
of  the  FORTRAN  statements  on  coding  sheets,  using  a  specific 


example. 
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SECTION  l 

PROBLEM  FORMULATION  AND  GOVERNING  EQUATIONS 

In  this  section,  a  typical  problem  in  subliving  ablation  will  be  formulated 
and  discussed  in  order  to  bring  out  the  sngineering  techniques  adopted.  In  the 
section  following  this  one  the  structure  and  processes  of  the  computer  program 
itself  will  be  described.  The  present  discussion  of  engineering  techniques 
culminates  in  the  formulation  of  a  forward-difference  numerical  procedure.  Along 
the  way,  various  inputs  to  this  numerical  procedure  are  introduced.  These  inputs 
are  reviewed  in  Section  2. 

In  general,  the  physical  problem  consists  of  a  flight  body  whose  surface  is 
subjected  to  aerodynamic  heating  while  also  exchanging  thermal  radiation  with  the 
surroundings  (i.e.,  "effective  space").  Figure  1  shows  schematic  diagrams  for  a 
typical  problem. 


✓ 


(a)  {b) 

Figure  1.  Schematic  Diagram  of  a  Typical  Ablation  Problem 

The  center  of  interest  in  the  overall  model  in  Figure  la  is  shown  as  a  segment 
of  the  boundary  layer  and  wall  in  Figure  lb.  In  order  to  determine  local  flow 
conditions  (Figure  lb),  the  aerodynamic  problem  (Figure  la)  must  first  be  solved. 
Completion  of  the  aerodynamic  phr.se  of  the  calculations  makes  it  possible  to  carry 
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out  the  second  phase  which  involves  the  determination  of  local  heat  and  mass 
transfer  at  the  exposed  surface.  The  second  phase  of  the  calculations  can  be 
carried  out  if  the  "present-time"  local  surface  temperature  has  b»rn  calculated 
or  the  "present-time"  ablation  surface  temperature  has  been  specified  on  the 
basis  of  published  ablation  data.  The  third  phase  of  the  problem  solution  makes 
use  of  finite-difference  equations  and  carries  heat  conduction  in  the  wall  from 
"present  time"  to'*future  time"  by  taking  a  forward  step  in  time. 

The  overall  ablation  problem  is  defined  by  setting  forth  body  configuration, 
initial  conditions,  boundary  conditions,  and  material  properties.  Problem 
definition  for  the  region  external  to  the  wall  surface  involves  the  following 
items : 

(1)  Flight  body  configuration  and  local  flow  conditions 

(2)  Mach  number  and  altitude  versus  time  after  launch 

(3)  Atmospheric  properties  versus  altitude 

(4)  Transition  Reynolds  number 

(5)  Radiation  environment 

Problem  definition  for  the  region  at  and  beneath  the  wall  surface  involves  the 
following  items. 

(1)  Wall  materials,  their  arrangement  and  dimensions 

(2)  Wall  material  properties  as  functions  of  temperature: 

a.  Thermal  conductivity 

b.  Density 

c.  Specific  heat 

d.  Surface  enissivity 

(3)  Additional  properties  for  ablators  as  functions  of  (15  -  H  ): 

r  aw  w 

a.  The  modified  effective  heat  of  afation 


b.  Ablation  surface  temperature 
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(4)  Initial  wall  temperature  distribution 

(5)  Thermal  boundary  conditions  on  internal  surfaces  of  the  vehicle 
Under  these  conditions  it  can  be  shown  that  the  wall  surface  temperature  and 
ablation  rate  are  determined  at  every  point  in  flight. 

In  the  example  of  the  flat  plate  at  angle-of-attack  (Figure  la),  computations 
begin  with  the  conversion  of  upstream  flow  conditions  to  the  local  conditions  of 
the  stream  external  to  the  boundary  layer  at  the  body  station  under  study.  Com¬ 
pressible  flow  charts  or  air  cables  are  used  to  determine  pressure  and  Mach  number 
behind  the  shock  using  the  perfect  gas,  the  thermally  perfect  gas,  or  the  real 
gas  case  ss  appropriate.  Next,  the  energy  equation  is  used  to  determine  the  local 
temperature.  At  this  stage  the  local  pressure,  temperature,  and  Mach  number  ac 
the  edge  of  the  boundary  layer  are  sufficient  to  determine  the  local  heat  transfer 
if  the  present-time  local  wall  surface  temperature  is  known.  Fortunately,  the 
latter  is  always  known  in  a  forward-difference  numerical  method  of  solution. 

This  completes  the  determination  of  local  flow  conditions.  One  additional  cal¬ 
culation  is  needed  to  determine  the  local  Reynolds  number. 

The  local  convective  heat  transfer  is  based  on  standard  equations  of  the  form 


he  x 

-V-  -  cV 
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C,  m,  and  n  are  determined  by  the  prescribed  geometry  and  by  the  relative  mag¬ 
nitudes  of  the  previously  calculated  local  Reynolds  number  and  the  prescribed 
transition  Reynolds  number.  During  ablation,  the  effect  of  mass  transfer  on 
convection  heat  transfer  must  be  taken  into  account.  This  effect  i  provided  for 
later  on  in  the  formulation  of  the  energy  equation  for  an  ablating  wall  surface 
element  by  introducing  the  transpiration  factor  3,..  Fcr  the  present,  the 

/I 

enthalpy-baeed  heat  transfer  coefficient  is  employed  in  rhe  relation 
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q  -  hu  (H  -  H  ) 
no  H,ov  aw  w 


(2) 


to  obtain  the  ordinary  local  rate  of  convective  heat  flow  per  unit  .area  to  a 
non-ablating  wall. 

The  convective  calculation  procedure  cills  for  determining  air  properties 
by  a  modification  of  Eckert's  (1)  reference  enthalpy  (i.e.,  K*)  method.  The 
modification  involves  the  uso  of  numerical  constants  recommended  by  Rubesin 
and  Johnson  (2)  and  Sotaner  and  Short  (3),  respectively,  as  follows 

H*  =■  0.58  H  +  0.42  (1  +  0.076  M2 )H  acinar)  (3) 


H* 


0,45  H  +  0.55  (1  +  0.054  M2)}!  (Turbulent) 


w 


(4) 


A  value  of  T*  corresponding  to  H*  is  found  from  a  polynomial  based  on  well- 
accepted  air  tables  of  H  versus  T  at  a  constant  pressure  of  one  atmosphere.  With 
T*  so  determined,  an  enthalpy-based  heat-transfer  coefficient  is  computed  u^ing 
a  rearranged  form  of  the  standard  heat  transfer  relation  containing  Pr  “  0.72. 

The  relation  is 


h,  -  ex""1  (p*  „>■  C.V  (0.J2)”-1  (5) 


The  Sutherland  viscosity  relation 


.*3/2 


O  ie 

T  +  X 


(6) 


provides  the  viscosity. 

The  adiabatic-wall  enthalpy  of  the  bemdary- layer  air  is  computed  by 
using  a  recovery  factor  of  1,0,  0.85,  or  0.9  for  the  stagnation  point,  the  .laminar 
boundary  layer,  or  the  turbulent  boundary  layer  respectively.  Local  convecliv^ 
(and  non-ablative)  heat  transfer  qQ  is  thin  computed  by  p-se  of  Equation  (2)  to 
complete  the  convective  phase  of  the  calculation. 
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Radiation  from  the  wall  surface  to  the  prescribed  radiation  environment  is 
computed  by  means  of  the  standard  radiation  equation  based  on  a  <  nit  area  of 
surface 


‘‘rad 


a  F.F  {T4 
A  e  w 


T4  ) 
space 


<7) 


Configuration,  emissivity  factors,  and  effective  space  temperature  Tgpace  used 
in  this  equation  are  prescribed  during  the  problem  definition  stage  and  the 
present-time  wall  temperature  is  always  known  in  the  forward-difference  method. 

The  variable-property  conduction  situation  in  the  wall  is  handled  by 
employing  a  lumped -parameter  equation  of  non-steady  heat  flow  in  more  than  one 
dimension.  A  two-dimensional  sketch  of  a  typical  lumped -parameter  break-up  of  a 
wall  segment  is  shown  in  Figure  2.  This  view  of  a  wall  segment  is  an  enlargement 
of  the  segment  shown  earlier  in  Figure  lb. 
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Figure  2.  A  Two-Dimensional  Break-up  of  a  Wall  Segment  Showing  Nodes 

In  the  absence  of  ablation,  convection  heat  transfer  per  unit  area  qQ  arrives 
at  area  A^  of  element  1  and  ■tra(jA^  departs.  Additional  heat  is  transferred  to 
element  1  by  conduction  from  other  elements  i  through  conductances  kA/f ,  k.^/L, 
etc.  which  are  called  K.  ..  The  net  heat  transfer  to  element  1  is  equal  to 
the  thermal  storage  in  a  thermal  capacitance  (pVc^,^),  »  C^.  Thus,  if  tempera¬ 
tures  of  elements  i  at  the  present  time  t  are  denoted  by  T^  and  those  at  future 
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time  t  1  it  are  denoted  by  1^,  a  forward-difference  equation  can  be  written  as 


<UA 


o"l 


E 


Ki,i^Ti  '  V  qradAl 


(8) 


To  prevent  temperature  overshoot  and  subsequent  oscillations  in  calculated 
values,  Dusinberre  (4)  shows  that  the  prescribed  time  step  At  must  be  less  than 
a  critical  time  step 


V  5*1.1 


(9) 


where  the  K.  .  include  the  respective  equivalent  conductances  for  convection 
i)* 

and  radiation, 


Vi 


T  ,  -  T. 
aw,l  1 


(10) 


and 
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.  time  step  siz~d  by  Dusinberre's  criterion  is  called  a  stability  time  step. 

When  an  element  is  not  subjected  to  convection  or  radiation  the  corresponding 

terms  drop  out  of  Bquations  (8)  and  (9). 

As  the  calculated  temperature  of  element  number  1  (Figure  2)  moves  toward 

the  ablation  temperature  of  element  number  1,  overshooting  and  oscillation  must 

be  avoided.  It  is  assumed  that  the  ablation  temperature  associated  with  the 

prescribed  material  of  the  protective  layer  has  been  specified  as  a  function  of 

heat-transfer  driving  potential  (H  -  H  )  or  as  a  constant.  Thus,  ablation 

aw  w 

chemistry  is  accounted  for  in  the  calculations  empirics. ly,  and  the  need  for 
specific  temperature-dependent  chemical  relations  is  circumvented.  With  ablation 
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temperature  specified  as  a  constant  or  in  terms  of  (H  -  H  ),  it  is  known  as 

aw  w 

the  calculations  move  ahead  in  flight  time.  Any  overshooting  of  the  ablation 
temperature  T^,  ^  by  surface  temperature  is  prevented  by  making  the  time 
step  in  equation  (S)  less  than  or  equal  to 


-  V 


fVi +  E  V(V  Ii)  •  ,.»dAi] 

i=2 


(12) 


In  the  calculation  procedure  this  criterion  is  tested  along  with  the  stability 
expression  (9)  and  the  smallest  of  the  two  time  steps  is  selected.  In  practice 
all  elements  are  scanned,  and  the  time  step  selected  is  the  smallest  for  all 
elements . 

If  surface  element  1  is  ablating,  the  ablation  processes  can  be  likened  to 
sublimation.  Accordingly,  the  ordinary  convective  heat  transfer  toward  area 
A^  is  diminished  by  the  equivalent  blocking  action,  ™Pjj(Kaw  j -  ^)  due  to  the 

mass  blowing  rate  m,  and  transpiration  factor  (see  Appendix  A).  Hence  qQ 
in  Equa'  ion  (8)  is  replaced  by 


qo  - 


&„(H  .  -  H  .) 

H  aw,l  w.l 


(13) 


during  ablation.  The  temperature-dependent  ablation  chemistry  is  represented, 
in  the  sublimation  model,  by  an  enthalpy  of  sublimation  and  a  sublimation 

.  temperature  T  .  If  ablation  occurs  during  a  forward  time  step,  the  ablating 
element  is  held  at  its  present  ablation  temperature  T  ^  until  the  completion 
of  the  time  step.  Therefore,  during  ablation  the  capacitance  term  of  Equation 
(8)  is  replaced  by  a  term  which  accounts  for  the  mas*'  decrease 


H(sg) 


P 


1 


(14) 
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Since  in  finite-difference  form  m  =  0,(Vj-  V^)/AjAt,  Equation  (8)  becomes,  upon 
introduction  or  Equations  (13)  and  (1A)  and  rearrangement  of  terms, 


Vi -  Z  Ki,i<v  Ti>  ■  ’„dAi 


H .  .  +  8U(H  -  H  ) 

(sg)  H  aw  w' 
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(15) 


Ihe  term  in  square  brackets  is  called  the  modified  effective  heat  of  ablation. 
In  Appendix  A  it  is  shown  that 


H,  ,  +  g„(H  -  H  )  = 

(sg)w  H  aw  w 


q  pt  q  J 
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semi- 
infinite 
test  body 


It  is  clear  thac  there  are  at  least  two  ways  of  obtaining  the  input  data  needed 
to  fill  in  the  square  bracket  in  Equation  (15).  They  are: 

(1)  Obtain  H,  .  from  the  chemical  literature  and  g„  from  the  mass-transfer 

(sg)  H 

literature,  and  evaluate  the  term  directly. 

(2)  Obtain  from  the  literature  on  ablation  experiments  thp  quantities 

needed  to  complete  the  semi-infinite  test  body  expression  in  Equation  (16). 

During  abla.ion  calculations ,  the  time  step  chosen  for  Equation  (15)  must 
not  be  such  that  becomes  negative,  for  this  would  mean  more  element  volume 
was  removed  than  was  available  at  the  beginning  of  the  time  step.  It  follows 
that  the  time-step  must  be  less  than  or  equal  to 


Qfai. :  y,.-  h„»i  <’ivi 
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(17) 
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In  the  final  analysis,  a  time  step  is  chosen  after  scanning  all  ablating  surface 
elements  to  check  for  the  smallest  time  to  ablate  away,  after . scanning  all  non¬ 
ablating  surface  elements  to  check  for  the  smallest  time  to  reach  ablation 
temperature ,  and  after  scanning  all  elements  to  check  for  the  smallest  stabill ty 
time  step.  The  smallest  time  step  dictated  by  all  three  criteria  must  not  be 
exceeded  by  a  given  step  forward  in  time. 

Figure  3  shoes  an  intermediate  stage  of  ablation.  Elements  are  assumed  to 
be  proportioned  so  that  ablation  of  a  given  element  occurs  in  one  dimension. 

As  a  consequence  of  the  effect  of  volume  change  on  the  cross-sectional  area  and 
path- length  of  conduction  paths  between  nodes,  adjustments  must  be  made  in  the 
K-values  of  Equations  (8),  (9),  (12),  (15),  and  (17).  With  reference  to  Figure  3 
the  thermal  conductance  between  nodes  i  and  1  is  defined  by  tne  relation 


i,l 


— i -  +  —  +  — 5— 

Ki,il  Kil  Kil,l 


(18) 


where  the  subscript  il  denotes  the  real  or  imaginary  interface  between  elements 
x  ana  1.  A  typical  conductance  term  is  defined  by  the  relation 
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i.il 


(19) 


where  a^^  is  the  area  normal  to  the  flow  or  heat  in  the  x-direction.  According 

d- 

to  Figure  3,  a^  must  be  continually  reduced  by  a  factor  V^/V^  in  order  to  account 

for  ablation.  Also,  l ^  ^  must  be  continually  reduced  by  the  same  ratio  to 

account  for  the  movement  of  node  i  toward  interface  i3.  By  means  of  these 

adjustments,  the  K,  .  in  Equations  (8),  (9),  (15),  and  (17)  are  modified  to 
1  j  X 

reflect  the  current  status  of  conduction  geometry. 
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Figure  3.  Surface  Element  Shrinkage  Due  to  Ablation 

With  regard  to  the  earlier  determination  of  the  local  convective  heat 
transfer  coefficient  h„  it  can  now  be  seen  that  by  taking  a  forward  step  in 
t ime,  the  calculation  always  has  available  for  use  the  local  wall  surface 
temperature  at  the  present  time.  Thus,  the  reference  temperature  T*  (see 
Equations  (3)  and  (4))  is  computed  for  conditions  at  the  present  time  and  in  the  f 

same  way,  material  properties  are  read  from  tables  of  properties  by  entering 
the  tables  at  point  temperatures  corresponding  to  the  present  time.  Lastly,  in  ■ 

the  forward  difference  method,  the  rate  of  heat  flow  to  any  node  (see  Equations  „ 

i 

(8)  and  (15))  is  always  based  on  temperatures  at  the  present  time.  The  method  A 

therefore  avoids  the  iteration  involved  in  a  formulation  calling  for  heat  ^ 

JL 

conduction,  convection,  radiation  or  property  inputs  corresponding  to  temperature 
conditions  at  the  end  of  the  time  step  under  consideration.  T 

I 

I 

1 

I 

i 
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SECTION  2 

STRUCTURE  OF  PROGRAM 

The  basic  ablation  computer  program  contains,  in  effect,  a  set  of 
available  thermal  capacitors  and  conductors  corresponding  to  capacitances  C 
and  conductances  K  of  the  previous  section.  These  are  connected  by  the  user 
who,  in  writing  a  jntrol  program,  brings  into  play  the  lumped-parameter 
equations  described  in  the  previous  section.  Convective,  conductive,  and/or 
radiative  heating  rates  are  computed  for  each  capacitor  in  accordance  with 
these  equations.  Each  thermal  capacitance  is  computed  from  the  si.se  and  thermal 
properties  supplied  as  inputs.  Prior  to  the  onset  of  ablation,  the  program 
confutes  transient  temperatures,  based  on  the  current  net  heating  rate  and 
thermal  capacitance,  by  taking  a  forward  step  in  time.  As  already  pointed  out, 
the  program  sizes  the  time  step  so  as  not  to  allow  uhe  temperature  of  any 
capacitor  to  overshoot  and  cause  oscillations  in  s^^sequent  finite  difference 
calculations  and  so  as  not  to  allow  the  temperature  of  any  surface  element  to 
exceed  its  ablation  teaqserature .  The  ablation  temperature  for  each  material 
is  an  input  specified  as  a  function  of  the  difference  bet*;een  the  adiabatic 
wall  and  wall  enthalpies  of  air.  When  an  element  reaches  its  ablation  temperature 
its  temperature  is  thereafter  determined  by  the  local  difference  between  the 
adiabatic  wall  and  wall  enthalpies  of  air.  If  the  heat  input  to  ar.  ablating 
element  becomes  negative  the  element  temperature  is  again  computed  from  the 
relation  between  net  heat  input  and  capacitance.  While  ablation  is  in  progress 
the  time  step  is  sized  to  prevent  the  calculated  volume  decrease  of  the  element 
from  overshooting  the  element  volume  available.  Also  provision  is  made  to 
vary  the  resistances  between  '-apacitors  (i.e.,  elements in  accordance  with  the 
decrease  in  surface  capacitor  volumes  due  to  ablation  and  to  discard  capacitors  which 
have  ablated  away.  ontained  in  the  printout  of  the  program  are  the  free  stream 
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aerodynamic  conditions,  the  local  flow  conditions  including  the  local  heat 
transfer  coefficient  and  adiabatic  wall  temperature,  the  capacitor  temperatures 
and  temperatures  at  the  interface  between  the  capacitors,  the  thickness  of 
remaining  ablator,  the  modified  effective  heat  of  ablation,  and  the  ablation 
temperature. 

The  ablation  program  has  been  constructed  in  two  distinct  sections;  a 
group  of  binary  subroutines  and  a  FORTRAN  control  program.  The  binary  subroutine 
deck  contains,  in  effect,  the  reservoir  of  capacitors  and  resistors  mentioned 
earlier.  It  also  contains  the  lumped -parameter  equations  necessary  to  compute 
the  various  modes  of  heat  transfer,  the  equations  to  compute  capacitor  tempera¬ 
tures,  and  equations  which  give  it  the  ability  to  size  the  time  step  according 
to  the  three  time-step  criteria.  The  geometry  being  analyzed  is  taken  into 
account  by  the  user  who  will  be  called  the  F0RTRAN  control  program  writer.  The 
F0RTRAN  control  program  writer,  in  affect,  assembles  a  network  from 
the  reservoir  of  capacitors  and  conductors  by  designating  specific  index  numbers 
in  the  calling  sequence  of  the  appropriate  binary  subroutines.  By  calling  other 
routines  he  se's  initial  conditions  and  provides  for  computation  of  heating  rates 
as  functions  of  time,  computation  of  transient  temperatures  and  ablation  rates, 
and  computation  of  properly  sized  time  steps.  This  feature  of  being  able  to 
construct  a  problem  in  the  course  of  calling  various  subroutines  gives  the  pro¬ 
gram  a  great  deal  of  flexibility. 

The  F0RTRAN  control  program  supplies  the  input  data  needed  by  the  subroutines 
to  make  the  required  calculations.  These  input  data  are  entered  into  the  orogram 
via  a  O0MM0N  statement  or  thiough  the  calling  sequence  (arguments)  of  the  sub¬ 
routine  CALL  statements.  Furthermore,  a  detailed  break-down  of  certain  of  these 
input  data  is  given  in  tabular  statementc. 
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Specifically,  the  F0RTRAN  control  program  consists  of  the  following  parts: 

I.  A  DIMENSION  statement 

II.  A  O0MM0N  statement 

III.  A  series  of  tabular  statements 

IV.  A  series  of  CALL  statements 

As  usual  the  Part  I  DIMENS10N  statement  sets  aside  the  appropriate  number  of 
tabular  spaces  in  the  computer  and  the  Part  II  O0MM0N  statement  provides  the 
appropriate  number  of  spaces  through  which,  in  effect,  certain  variables  can 
pass  from  one  subroutine  to  another.  A  small,  but  not  negligible,  amount  of 
information  on  problem  geometry  enters  the  O0MM0N  statement.  By  contrast,  the 
tabular  statements  of  Part  III  deal  almost  wholly  with  details  of  the  problem 
at  hand.  Tabular  statements  must  include  the  following: 

1.  Mach  number  versus  time 

2.  Altitude  versus  time 

3.  Ambient  temperature  and  pressure  versus  altitude 

4.  Local  aerodynamic  flow  conditions  (M/Mq  versus  Mo  and  P/Pq 
\ ersus  Mo) 

5.  The  effective  space  temperature  versus  altitude 

6.  The  emissivity  of  external  surfaces  versus  temperature 

7.  A  series  of  coefficients  and  exponents  concerned  with  the 
aerodynamic  heating  equations,  plus  a  transition  Reynolds  number 

8.  The  initial  volumes  of  all  potential  ablating  capacitors 

9.  An  ablation  sequence  prescribing  the  order  of  capacitor  ablation 

10.  The  material  density,  specific  heat,  aid  conductivity  as  a 


function  of  temperature 
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11.  The  ablation  temperature  versus  air-enthalpy  difference 

12.  The  modified  effective  heat  of  ablation  versus  air-enthalpy 
difference 

13.  The  capacitor  numbers  of  all  initial  surface  elements 

In  supplying  the  foregoing  tabular  information  in  numerical  form  the  F0RTRAN 
control  program  writer  must  assign  a  variable  name  to  each  type  of  entry  and 
identify,  by  subscripting,  each  individual  numerical  entry.  With  the  exception 
of  items  8  and  9,  Part  III  tabular  information  is  transferred  to  subroutines  via 
the  Part  IV  calling  sequences  (arguments)  of  the  subroutine  CALL  statements. 

This  transfer  is  made  possible  when  the  program  writer  enters  tabular  variable 
names  in  designated  positions  within  the  argument  cf  the  CALL  statements.  Items 
8  and  9  are  entered  into  the  O0MM0N  statement  and  thereby  become  directly  avail¬ 
able  to  the  subroutines.  In  writing  the  CALL  statements  of  Part  IV,  the  F0RTRAN 
control  program  writer  draws  on  a  reservoir  of  capacitors  and  conductors  and 
assembles  them  into  a  network.  This  step  is  describe-’  in  more  detail  in  the  next 
paragraph. 

Problem  inputs  not  covered  by  items  1  to  13  above,  are  introduced  into  the 
CALL  statement  arguments  along  with  the  tabular  variables.  They  appear  in  the 
argument  as  numerical  values,  unlike  the  tabular  variables  which  appear  as  varia¬ 
ble  names.  The  numerical  inputs  in  question  are; 

1.  Index  numbers  of  capacitors 

2.  Initial  capacitor  temperatures 

3.  Exposed  surface  area  of  capacitor  and  its  Reynolds  number  refer¬ 
ence  length 

4.  Capacitor  volume 
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5.  Original  area-to-length  ratio  for  resistors  and  contact  resistance 
’•alues 

6.  Radiation  configuration  factors 

It  is  important  to  note  again  that  the  flexibility  of  the  program  exists  in  the 
reservoir  of  capacitors  and  conductors  whose  index  numbers  can  be  called  according 
to  the  requirements  of  the  geometry  of  the  problem  under  study.  Experience 
teaches  one  how  best  to  divide  a  given  heat-transfer  model  into  the  elements 
that  ultimately  become  the  capacitors  and  conductors  of  the  computer  program. 

Under  Part  IV  of  the  F0RTRAN  control  program  twelve  main  subroutines  in  the 
binary  subroutine  deck  can  be  called.  They  are  entitled: 


1. 

ASET 

7. 

C0N 

2. 

TRAJ 

8. 

AO0MCN 

3. 

AIM 

9. 

RAD 

4. 

F0RALT 

10. 

WRITE 

5. 

AAER0 

11. 

ASTEP 

6. 

ACAP 

12. 

ABL 

For  some  of  these  subroutines,  the  order  in  which  the  writer  of  the  F0RTRAN 
control  program  calls  them  is  critical;  for  others  it  is  not.  To  avoid  errors, 
it  is  suggested  chat  the  subroutines  be  called  in  the  order  given  above  in  all 
cases.  Even  when  one  is  well  versed  in  r'ne  use  of  the  program  he  will  find  no 
particular  advantage  in  changing  the  order  of  the  CALL  statements  from  that 
given  above.  A  discussion  follows  which  is  intended  to  provide  a  brief  view  of 
the  structure  and  operation  of  each  s>  ’-routine  and  the  overall  operation  of  the 
program. 

The  ASET  subroutine  is  a  FAP  subroutine  which  sett  the  initial  temperature 
of  each  capacitor  and  defines  the  beginning  and  end  of  f:ight  time.  These 
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settings  are  chosen  by  the  F0RTRAN  control  program  writer  and  entered  into  the 
argument  of  t!:i  CALL  ASET  statement. 

Subroutine  TRAJ  is  essentially  a  table  look*up  F0RT5AN  routine  which 
requires  as  inputs,  through  its  argument,  tlic  variable  names  previously  assigned 
by  the  control  program  writer  in  preparing  the  tables  of  Mach  number  and  altitude 
versus  time.  The  function  Of  XRAJ  is  to  quadratically  interpolate  in. the  fables 
using  the  current  value  of  time  to  obtain  the  current  Mach  number  and  altitude. 
Tiese  two  items  leave  the  subroutine  via  the  C0MM0N  statement  and  thus  become 
available  for  later  use  by  the  ATM,  AAER0,  F0RALT,  and  WRITE  subroutines. 

The  F0RTRAN  subroutine  AIM  is  also  a  table  look-up  routine.  The  control 
program  writer  must  place  the  variable  names  assigned  in  the  tables  of  static 
temperature,  pressure,  and  altitude  in’ the  argument  of  ATM  when  calling 
the  routine.  ATM  linearly  interpolates  in  the  tables  using  the  current  altitude, 
supplied  by  TRAJ,  to  obtain  the  current  atmospheric  temperature  and  pressure  Pq. 
Both  items  are  then  supplied  to  AAEK0  and  WRITE  through  O0MM0N. 

F0RALT  is  a  F0RTRAN  subroutine  whose  function  is  to  quadratically  interpolate 
in  the  table  of  effective  space  temperature  versus  altitude  (using  the  altitude 
supplied  by  TRAJ)  tc  obtain  the  current  space  temperature.  The  routine  then 
assigns  this  temperature  to  a  particular  space  node.  When  calling  the  routine, 
the  control  program  writer  must  include  in  the  argument  of  the  CALL  statement, 
the  node  index  number  he  desires  to  represent  space  and  the  previously  assigned 
variable  names  for  the  effective  space  temperature  and  altitude.  The  space 
temperature  is  supplied  to  RAD  and  WRITE  through  a  O0MM0N  statement. 

The  AAER0  subroutine  is  a  FAP  subroutine.  Its  primary  function  is  to 
compute  the  aerodynamic  heating.  In  order  to  accompli,  this,  the  AAER0  FAP 
routine  is  broken  into  two  F0RIRAN  subroutines,  A£R0A  and  AAER0B.  The  control 


TSa  Mk;  Ui»l*«f».(y 

nnw  unuw> 

Wvw  feHn*,  Mvylan4 


21. 


program  writer  calls  AAER0  which  then  automatically  calls  AER0A  and  AAER0B. 

The  AER0A  portion  of  the  AAER0  routine  computes  an  adiabatic  wall  temperature 
and  an  "effective"  heat  transfer  coefficient  which  excludes  a  major  length 
dimension  termed  the  reference  length.  Information  needed  by  AER0A  is  supplied 
either  through  a  O0MM0N  statement  from  other  subroutines  or  by  the  control  pro¬ 
gram  writer  by  properly  satisfying  the  argument  when  calling  the  routine.  The 
control  program  writer  must  enter  in  the  calling  sequence  of  AAER0  the  variable 
names  previously  assigned  in  preparation  of: 

a.  The  table  giving  the  transition  Reynolds  number  plus  coefficients 
and  exponents  needed  in  the  aerodynamic  heating  equations 

b.  The  two  tables  giving  the  local  flow  conditions,  M/Mq,  etc. 

c.  The  two  tables  giving  the  modified  effective  heat  of  ablation 
and  ablation  temperature  versus  enthalpy  difference 

Also  there  must  be  included  a  series  of  values  for  the  index  number,  surface  area, 

and  reference  length  of  all  those  surface  capacitors  subjected  to  the  aerodynamic 

flow  conditions  listed  in  the  calling  sequence. 

Briefly,  he  AER0A  subroutine  operates  in  the  following  manner.  The  ratios 

of  local  Mach  number  to  free  stream  Mach  number  (M/M^)  and  local  pressure  to 

free  stream  pressure  (P/Pq)  as  a  function  of  free  stream  Mach  number  are  obtained 

through  the  argument  from  the  list  of  tables  in  the  control  program.  The  free 

stream  Mach  number  Mq  computed  by  TRAJ,  is  obtained  from  O0M0N  and  a  quadratic 

interpolation  is  used  to  obtain  the  current  values  of  M/M  and  P/P  .  Hie  local 

o  o 

Mach  number  and  pressure  are  obtained  by  multiplying  these  ratios  by  Mq  aid  P^ 
respectively,  where  the  value  of  Pq  is  obtained  from  ATM  through  O0JH0N.  Next, 
flight  and  local  Mach  numbers  are  used  in  an  energy  equation  to  convert  the 
temperature  ahead  of  the  vehicle  to  the  local  temperature  of  the  flow  external  to 
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the  boundary  layer.  A  local  Reynolds  number  is  then  computed  based  on  the  local 
temperature,  local  pressure,  local  Mach  number,  and  the  first  ’-eference  length 
given  in  the  AAER0  calling  sequence.  This  Reynolds  number  is  compared  with  a 
transition  Reynolds  number  supplied  in  the  list  of  tables  in  the  control  program 
to  determine  whether  subsequent  calculations  should  be  for  a  laminar  or  turbulent 
boundary  layer.  This  Reynolds  number  is  not  used,  however,  to  compute  the  heat 
transfer  coefficient.-  Instead,  a  reference  enthalpy  and  adiabatic  wall  enthalpy 
are  computed.  These  enthalpies  are  converted  to  temperatures  by  means  of  a 
polynomial  which  has  been  fitted  to  conventional  temperature-enthalpy  data  for 
air  at  standard  atmospheric  pressure.  The  adiabatic  wall  temperature  is  placed 
in  O0MM0N.  Air  properties  are  then  evaluated  at  the  reference  temperature,  and 
subsequently  combined  with  the  local  pressure,  local  velocity,  and  coefficients 
and  exponents  given  in  the  table  listings  in  the  control  program,  to  obtain  an 
"effective"  local  heat  transfer  coefficient.  Again,  the  word  "effective"  is  used 
to  indicate  that  a  reference  length  has  not  yet  been  included  in  the  computations. 
This  completes  the  AER0A  calculations  and  the  program  automatically  proceeds  to 
.vAER0B. 

The  primary  function  of  AAER0B  is  to  compute  the  aerodynamic  heating  rate  to 
each  capacitor  listed  in  the  calling  sequence  of  AABR0.  In  order  to  do  this  it 
applies  the  effective  local  heat  transfer  coefficient  to  each  capacitor  after 
obtaining  reference  lengths  from  the  AAER0  calling  sequence,  thereby  computing  a 
different  local  heat  transfer  coefficient  for  each  capacitor.  Next,  the  capacitor 
"surface  area"  is  obtained  from  tl  2  AAER0  calling  sequence,  the  local  adiaKstic 
wall  temperature  is  obtained  from  AER0A  through  O0MM0N,  and  the  capacitor  "surface 
temperature"  is  obtained  from  ASTEP  through  C0MM0N  frori  the  previous  time.  AAER0B 
then  has  all  the  information  needed  to  compute  the  heating  rate  to  each  capacitor. 
Tills  heating  rate  is  then  placed  in  O0MM0N. 
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A  secondary  function  of  AAER0B  is  to  determine  the  current  value  of  the 

ablation  temperature,  TABL,  and  modified  effective  heat  of  ablation,  HEFF.  In 

order  to  do  this,  it  obtains  the  current  value  of  the  enthalpy  difference, 

H  -  H  from  AER0A  and  uses  this  current  value  to  quadratieally  interpolate  in 
aw  v  *“ 

the  tables  of  HEFF  and  TABL  versus  H  -  H  .  As  previously  mentioned,  the  vari- 

aw  w 

able  names  assigned  to  these  tables  must  be  written  into  the  AAER0  calling  sequence. 

The  function  of  the  F0RTRAN  subroutine  ACAP  is  to  compute  the  thermal  capa¬ 
citance  of  each  capacitor.  Tabular  values  of  material  density  and  specific  heat 
are  given  in  the  control  program  as  functions  of  temperature  and  are  entered  into 
the  routine  by  listing  the  variable  names  of  the  tables  in  the  calling  sequence. 
Also,  the  capacitor  volume  is  entered.  The  routine  determines  the  current  density 
and  specific  heat  by  quadratic  interpolation  in  the  tables,  using  the  capacitor 
temperature  obtained  through  O0MMON  from  the  previous  time.  The  density  and 
specific  heat  are  multiplied  by  the  capacitor  volume  to  obtain  the  thermal  capa¬ 
citance  which  is  then  placed  in  O0MM&  One  call  of  subroutine  ACAP  handles  only 
one  capacitor.  Therefore,  the  control  program  must  call  ACAP  as  many  times  as 
there  are  number  of  capacitors.  It  is  worth  noting  that  by  repeatedly  calling 
index  numbers  in  the  ACAP  subroutine,  the  F0RT*iAN  control  program  writer  draws 
on  a  reservoir  of  capacitors  and,  in  essence,  defines  a  set  cf  nodes. 

F0RTRAN  subroutine  O0N  is  one  of  the  two  routines  in  the  pri  jram  which  compute' 
the  thermal  conductance  between  nodes  and  the  subsequent  conductive  heating  rate 
between  nodes.  By  introducing  appropriate  pairs  of  capacitor  index  numbers  into 
the  calling  sequence  of  the  O0N  routine  the  control  program  writer  draws  upon  a 
reservoir  of  conductors  and  assembles  a  node  network  into  a  conductor-capacitor 
configuration  approximating  the  thermal  model.  This  s  broutine  can  be  called 
repeatedly  to  join  together  two  nodes  of  non -ablating  materials.  When  O0N  is 
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called,  the  control  program  writer  lists  in  the  calling  sequence  the  two 
capacitor  numbers  being  joined,  the  ratio  of  area  to  length  between  nodes,  and 
the  variable  names  previously  assigned  in  preparing  the  tables  of  thermal  con¬ 
ductivity  versus  temperature.  The  routine  takes  this  information,  refers  to 
C0MM0N  for  the  node  temperatures  at  the  end  of  the  previous  time  step,  and 
computes  the  heating  rate.  This  heating  rate  is  then  added  algebraically  to 
the  heating  rate  (if  any)  previously  computed  by  AAER0.  At  this  stage,  there¬ 
fore,  the  heating  rate  stored  in  O0MM0N  includes  both  aerodynamic  and  conductive 
heating  rates  on  capacitors  designated  by  the  program  writer.  Optionally,  the 
program  writer  may  have  conduction  heat  transfer  computed  by  subroutine  AO0MCN 
as  described  in  the  next  paragraph. 

The  F0RTRAN  subroutine  AO0MCN  is  a  routine  which  computes  the  thermal 
resistance  between  ablative  as  well  as  non-ablative  capacitors.  To  obtain  the 
overall  thermal  resistance  between  two  nodes,  three  separate  resistances  are 
computed;  one  from  the  center  of  one  node  to  the  interface,  a  second  across  the 
interface,  and  a  third  from  the  interface  to  the  center  of  the  second  node.  This 
feature  makes  .  G0MCN  handy  to  use  in  cases  where  either  the  material  conduction 
area  changes  at  the  interface  or  it  is  desired  to  account  for  interface  resistance. 
AC0MCN  also  computes  temperatures  on  each  side  of  the  interface  and  feeds  this 
information  through  O0MM0N  to  the  WRITE  subroutine.  (It  should  be  noted  here 
that  all  other  temperatures  are  computed  by  the  ASTEP  subroutine.)  When  calling 
AO0MCN,  the  F0RTRAN  control  program  writer  must  list  the  following  items  in  the 
calling  sequence: 

a.  The  index  numbers  of  the  two  capacitors  being  joined 

b.  The  area  to  length  ratio  from  the  first  node  to  the  interface 

c.  The  area  of  the  interface 
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d.  The  area  to  length  ratio  from  the  interface  to  the  second  node 

e.  T..e  variable  name  previously  assigned  to  the  tables  of  thermal 
conductivity  versus  temperature  for  the  material  renresented  by 
the  first  node 

f.  The  variable  name  previously  assigned  to  the  fables  of  thermal 
conductance  versus  temperature  of  the  interface 

g.  The  variable  name  previously  assigned  to  the  tables  of  thermal 
conductivity  versus  temperature  for  the  material  represented  by 
the  second  node 

Prior  to  any  computations  of  heating  rate  or  interface  temperature,  AO0MCN  first 
examines  a  counter  (which  is  set  in  the  ABL  routine)  to  determine  whether  either 
of  the  two  capacitors  being  joined  is  ablating.  If  neither  capacitor  is  ablating, 
the  subroutine  computes  the  interface  temperature  and  heating  rate  to  both 
capacitors.  For  this  purpose  the  subroutine  obtains  the  area  to  length  ratios 
and  interface  area  directly  from  the  calling  sequence,  the  thermal  conductivities 
from  the  calling  sequence,  and  present  capacitor  temperatures  from  OpHnpii.  The 
interface  temperature  is  then  computed  and  placed  in  O0MM0N.  The  new  increment  of 
heating  rate  for  each  of  the  two  indices  being  joined  is  added  algebraically  to 
the  amount  already  in  O0MM0N.  However,  if  after  the  counter  is  examined ,  one  or 
both  of  the  capacitors  are  found  to  be  ablating,  either  an  area  or  length  between 
nodes  must  be  adjusted  to  account  for  the  effect  of  decreasing  element  volume. 
Since  each  ablation  sequence  in  the  control  program  lists  capacitor  node  index 
numbers  in  a  column  sequence  starting  at  the  surface  and  extending  into  the 
ablator,  the  sequence  is  examined  by  the  routine  to  see  if  both  capacitor  index 
numbers  appear  in  a  given  sequence.  If  they  do,  the  capacitors  must  lie 
above  one  another.  If  only  one  is  found  in  a  given  ablation  sequence,  it  follows 
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that  the  two  capacitors  are  beside  each  other.  Once  the  routine  has  determined 
the  relative  location  of  the  capacitors,  either  the  original  area  or  original 
length  is  decreased  by  a  factor  equal  to  the  ratio  of  current  capacitor  volume 
(obtained  from  the  ABL  routine  through  C0MM0N)  to  original  volume  (obtained  from 
a  table  in  the  control  program  through  O0MM0N).  Once  the  area  to  length  ratios 
have  been  adjusted,  the  subroutine  proceeds  in  a  fashion  identical  to  the  non¬ 
ablating  case. 

F0RTRAN  subroutine  RAD  computes  the  radiation  heat  transfer  rate  between  the 
capacitor  numbers  listed  in  its  calling  sequence.  The  roucine  obtains  the  sur¬ 
face  ?r„a  directly  from  the  calling  sequence  and  obtains  the  surface  emissivity 
from  the  emissivity- temperature  table  supplied  in  the  control  program  whose 
variable  name  is  also  listed  in  the  calling  sequence.  The  wall  temperature  from 
the  end  of  the  previous  time  step  is  used  to  quadratically  interpolate  in  the 
table  and  obtain  the  current  emissivity.  The  emissivity  and  wall  temperature  are 
combined  with  the  effective  space  temperature  obtained  from  F0RALT  through  C0MM0N, 
to  calculate  the  radiative  heating  rate  on  the  surface  element  in  question.  The 
net  heating  ra!  e  for  each  capacitor  surface  element  is  then  obtained  from  O0MM0N, 
altered  by  the  amount  of  the  radiative  heating  rate,  and  replaced  in  O0MM0N. 

The  sole  function  of  the  FAP  routine  WRITE  is  to  handle  the  printout  of  the 
data.  The  WRITE  subroutine  itself  calls  a  F0RTRAN  routine  entitled  PRNTA.  PR.NTA 
contains  a  predetermined  printout  format  which  is  executed  at  the  flight  times 
specified  by  the  control  program  writer  in  the  calling  sequence  of  WRITE.  At 
each  of  the  specified  flight  times,  WRITE  will  print  out  the  following  ite^s : 

1.  Flight  time,  altitude,  free  stream  Mach  number,  free  stream  pressure, 
and  free  stream  temperature 
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2.  Capacitor  temperatures  (if  they  have  changed  from  the  initial  values 
specified  in  the  argument  of  the  CALL  ASET  statement) 

3.  Contact  (interface)  temperatures 

k.  Local  flow  conditions;  Mach  number,  pressure,  temperature,  Reynolds 
number,  adiabatic  wall  temperature,  and  the  "effective"  heat  transfer 
coefficient 

5.  The  index  number  of  the  capacitor  nearest  to  its  ablation  temperature 

and  the  time  needed  for  it  to  reach  ablation  temperature 

6.  The  index  number  of  the  capacitor  nearest  to  being  completely  ablated 

away  and  the  time  to  complete  ablation 

7.  The  stability  time  step  and  the  capacitor  ntmber  associated  with  it 

8.  The  capacitor  number,  current  volume,  modified  effective  heat  of 
ablation,  and  ablation  temperature  for  all  capacitors  currently  ablating 

In  addition*To  printing  out  these  items  at  the  times  called  for  in  the  calling 
sequence,  the  program  will  print  out  any  time  a  capacitor  either  reaches  its 
ablation  temperature  or  completely  ablates  away,  when  ablation  is  interrupted  by 
a  negative  heating  rate,  or  when  ablation  is  resumed  following  a  return  to  a 
positive  heating  rate.  Such  printouts  are  triggered  by  a  code  supplied  by  the  ABL 
routine. 

The  function  of  F0RTRAN  subroutine  ASTEP  is  to  properly  size  the  time  step 
and  compute  new  capacitor  temperatures.  All  the  information  needed  by  ASTEP  to 
perform  calculations  is  obtained  from  other  routines.  It  has  no  calling 
sequence.  ASTEP  determines  three  different  ' ime  steps: 

(a).  The  miuimum  of  the  times  needed  by  each  surface  capacitor  to  reach 
its  ablation  temperature.  See  expression  (12) 
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(b) .  The  minimum  of  the  times  needed  by  each  surface  capacitor  to  reach 

complete  ablation.  See  expression  (17) 

(c) .  The  stability  time  step  needed  to  prevent  oscillation.?  in  subsequent 

finite  difference  calculations.  See  expressic  (9) 

ASTEP  chooses  the  minimum  of  (a),  (b),  and  (c)  to  step  ahead  in  time.  Time  step 
(a)  is  computed  using  the  capacitance  supplied  by  ACAP,  the  ablation  temperature 
supplied  by  AAER0B,  the  wall  temperature  from  the  previous  time,  and  the  not 
heating  rate  computed  by  AAER0,  AO0MCN,  C0N,  and  RAD.  Time  step  (b)  is  computed 
using  the  density  from  ACAP,  the  current  volume  from  ABL,  the  modified  effective 
heat  of  ablation  from  AAER0B,  and  the  net  heating  rate  from  AAER0,  AO0MCN,  C0N, 
and  RAD.  Time  step  (c)  is  computed  from  the  capacitance  supplied  by  ACAP  and  the 
cumulative  thermal  resistances  comouted  in  AAER0B,  O0N,  AO0MC*I,  and  RAD.  Follow¬ 
ing  selection  of  the  minimum  time  step,  the  routine  computes  new  temperatures 
based  on  the  previous  temperature  of  the  capacitor,  the  net  heating  rate,  the 
thermal  capacitance,  and  the  minimum  time  step. 

The  last  routine  is  ABL,  a  FaP  routine  which  calls  the  F0RTRAN  routine  ABL1. 
The  function  of  ABL1  is  to  compute  the  current  volume  of  ablating  capacitors  and 
discard  those  capacitors  which  have  completely  ablated  away.  It  also  sets  a 
counter  so  other  subroutines  will  know  whether  or  not  a  par  icular  capacitor  is 
currently  ablating  or  is  completely  ablated  away.  Finally,  j  changes  the  surface 
capacitor  when  an  element  completely  ablates  away  so  th..  convect../.i  (AAF.R0)  and 
radiation  (RAD)  will  be  applied  to  the  newly  exposed  capacitor.  When  calling 
the  subroutine,  the  writer  must  include  in  »-he  calling  sequence  the  variable 
names  assigned  to  the  numbers  of  the  initial  surface  capacitors.  These  variable 
names  and  associated  capacitor  numbers  must  have  bfce.i  {  -eviously  defined  in  the 
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control  program.  All  other  information 
The  future  volume  is  computed  using  the 
the  minimum  time  step,  the  density,  and 


needed  by  ABL  is  obtained  through  C0MM0N. 
present  volume,  the  net  heating  rate, 
the  modified  effective  hc.at  of  ablation. 
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SECTION  3 

THE  F0RTRAN  CONTROL  PROGRAM 

The  preceding  discussion  of  the  computer  program  has  laid  the  groundwork 
for  giving  specific  instructions  to  the  user.  This  section  identifies  the 
F0RTRAN  control  program  statements  and  gives  the  nomenclature.  As  noted  in 
the  previous  section,  the  F0RTRAN  control  program  consists  of: 

I.  DIMENSI0N  statements 

II.  A  C0MM0N  statement 

III.  A  series  01  tabular  statements 

IV.  A  series  of  CALL  statements 
These  statements  have  the  following  form  and  composition: 

I.  DIMENSI0N  Statements 

DIMENSI0N  F1(I),  F1V(J),  F2(K), . 

FI (I)  -  dependent  variable  table  Fl,  I  spaces  needed 

FlV(J)  -  independent  variable  table  F1V.  J  spares  needed 
F2(K)  -  dependent  variable  table  F2,  K  spaces  needed 

F2V(L)  *  Independent  variable  table  F2V,  L  spaces  needed 
Sl(K)  -  etc. 

SlV(K)  -  etc. 
etc. 


It  should  be  emphasized  that  these  variable  names  are  to  be  devised  by 
the  control  program  writer.  Floating  point  names  are  to  be  used. 
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DIMENSION  1ASEQ(500),  V(1000),  DUM<4407),  1AB(500)S  V0(1OOO),  DAMN<4002) 
lASoQCSOO'i  -  ablation  sequence  ,  500  spaces  available 
V(1000)  -  volume  of  capacitors,  1000  spaces  available 

DUM(4407)  "  44 07  dununy  spaces 

IAB(500)  -  ablation  "on-off"  index,  500  spaces  available 

V0(1OOO)  -  initial  volume  of  potential  ablation  elements,  1000 
spaces  available 
DAMN(4002)  -  4002  dummy  spaces 
II.  O0MH0N  Statements 

COMMON  T,  Q,  DUM,  IAB,  IASEQ,  DAMN 
T  -  capacitor  temperature 

Q  -  accumulated  heating  rate  being  applied  to  capacitor 
DUM  -  dummy 

1AE  -  defined  in  part  I 

IASEQ  -  defined  in  part  1 

V0  -  defined  in  part  1 

DAMN  -  dummy 

III.  Tabular  Statements 

Fl(l)  a  a  is  the  first  value  of  the  dependent  variable  in  a  table 
named  FI 

F1  (2)  "  b  b  ia  tbe  second  valua  of  the  dependent  variable  in  a  table 
named  FI 

PI(n)  »  k 

F1V(1)  *  n  n  is  the  total  number  of  items  in  ia''le  FI 
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F1V(2)  *  a'  a'  is  the  value  of  the  independent  variable  corresponding 
to  the  dependent  variable  Fl(l) 

It  should  be  noted  again  that  these  variable  names  are  to  be  devised  by  the 
control  program  writer. 

IV.  CALL  Statements 

1.  CALLASET (START,  ST0P ,  TEMPIN,  I,  INDEXI ,  Tl... INDEXI,  T(I)) 

START  -  time  at  beginning  of  computations 
ST0P  -  time  at  end  of  computations 

TEMPIN  -  temperature  at  which  certain  capacitors  are  to  be  initialized 
I  -  number  of  capacitors  not  desired  to  be  initialized  at 

TEMPIN 

INDEXI  -  capacitor  numbers  whose  initial  temperature  is  not  desired 
to  be  TEMPIN 

T(I)  -  initial  temperature  of  INDEXI 

2 .  CALLTRaJ  (F0FXM , XM , FOFXA , XA ) 

Ffl^XM  -  variable  name  given  to  free  stream  Mach  number  table 

XM  -  variable  name  given  to  corresponding  flight  time  table 

F0FXA  -  variable  name  given  to  altitude  table 

XA  -  variable  name  given  to  corresponding  flight  time  table 

(The  variable  names  are  actually  to  be  entered  in  the  form  deyised 
in  setting  up  each  table  in  part  III.) 
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3.  CALLATM  (F0FXT0,XT0,F0FXP0,XP0) 


4. 


r^Fxv0 


XT0 

F0FXP0 

XPC 


-  variable  name  given  tc  ambient  temperature  table 

-  variable  name  given  to  corresponding  altitude  table 
variable  name  given  to  ambient  pressure  table 

-  variable  name  given  to  corresponding  altitude  table 


CALLF0RALT  (INDEX,  F0FX,  X) 


INDEX  -  number  of  the  capacitor  designated  to  represent  space 


F0FX 


-  variable  name  given  to  rpace  temperature  table  (i.e.,  the 


table  gives  the  schedule  of  temperatures  to  be  assigned 
to  INDEX) 


X  -  variable  name  assigned  to  altitude  table  corresponding 
to  F0FX 

5.  CALLAAER0(ID,I,FL0N0S ,FIFXLM,XLM,F0FXLP,XLP,ABLT,HDIF1 ,HEFF ,KDIF2 , 

INDEX  1,  GE0N01,  P0SN01, . INDEXI,  GE0N0I,  P0SN0I) 

ID  aerodynamic  block  number  (a  block  being  a  group  of  surface 

caoacitors  all  imrior  t-V»£»  A.U  -  - -  1.  1 

“  “  -  - -  '•■**  V/J.  vuc  OCU1C  i-UCcti. 

flow  conditions) 

I  -  number  of  surface  capacitors  in  the  block 
FIX*N0S  -  variable  name  given  to  the  table  specifying  the  transition 
Reynolds  number  and  the  coefficients  and  exponents  needed 
to  compute  the  local  heat  transfer  coefficiei.t 
F0FXLM  -  variable  name  given  to  the  table  of  the  local-to-free-streero 
Mach  number  ratios 


XLM 


-  variable  name  given  to  corresponding  free  stream  Mach  number 


table 
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F0FXLP  -  variable  name  given  to  the  table  of  local-to-free-stream 
static  pressure  ratios 

XL?  -  variable  name  given  to  corresponding  free  stream  Mach 
number  table 

ABLT  -  variable  name  given  to  the  ablation  temperature  table 

HDIF1  -  variable  name  given  to  corresponding  air-enthalpy 
difference  table  (HAW  -  HW) 

HEFF  -  variable  name  given  to  the  modified  effective  heat  of 
ablation  table 

HDIF2  -  variable  name  given  to  corresponding  air-enthalpy  difference 
table 

INDEX1  -  capacitor  number  of  the  first  capacitor  in  the  block 
GK0N01  -  surface  area  of  first  capacitor 
P0SN01  -  reference  length  of  first  capacitor 
INDEXI  -  capacitor  number  of  the  last  capacitor  in  the  block 
GE0N0I  -  surface  area  of  last  capacitor 
P01N0I  -  reference  length  of  the  last  capacitor 
6.  CALLACAP  (INDEX.GE0N0,R,XR>CP,XCP) 

INDEX  -  capacitor  number 

GE0N0  -  capacitor  volume 

R  -  variable  name  given  to  the  table  listing  the  density 

of  the  capacitor  material 

XR  -  variable  name  given  to  the  corresponding  temperature  table 

CP  -  variable  name  given  to  the  table  listing  the  specific 

heat  of  the  capacitor  material 

XCP  -  variable  name  given  to  corre°pondir.g  temperature  table 
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7.  CALLC0N(INDEX1 ,  INDEX2 ,  GE0N0,  FOFX,  X) 

Ili'jEXl ,  INDEX2  -  identifying  numbers  of  the  two  capacitors 
being  joined 

GE0N0  -  the  trea  to  length  ratio  between  the  two  capacitors 
F0FX  -  the  variable  name  assigned  to  the  table  giving  the  thermal 
conductivity  of  the  material 

X  -  the  variable  name  given  to  the  corresponding  temperature  table 

8.  CALLAC0MCN(INDEX1 ,  INDEX2,  GE0N01,  GE0N02 .  GE0N03  ,  F0FX1 ,  Xl , 

F0FX2,  X2,  F0FX3 ,  X3) 

INDEXl ,  INDEX2  -  identifying  numbers  of  the  two  capacitors 
being  joined 

GE0N01  -  the  ratio  of  area  to  length  between  the  node  of  INDEXl 
and  the  interface 

GE0N02  -  the  area  of  the  interface 

GE0N03  -  the  ratio  of  area  to  length  between  the  interface  and  the 

r.CuC  VSJ.  XHUUA4 

F0FX1  -  variable  name  given  to  the  table  listing  the  thermal 
conductivity  of  the  material  represented  by  INDEXl 

XI  -  variable  name  of  the  table  listing  the  corresponding 

temperature 

F0FX2  -  the  variable  name  given  to  the  table  listing  the  interface 
conductance 

X2  -  variable  name  of  the  corresponding  temperature  table 

F0FX3  -  variable  name  given  to  the  table  listing  the  thermal 
conductivity  of  the  material  repit-'ented  by  INDEX2 
X3  -  variable  name  of  the  corresponding  temperature  table 
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9.  CaLLRAD(INLEX1,  INDEX2 ,  GE0N0,  F0FX,  X) 

INDEX. ,  INDEX2  -  identifying  numbers  of  the  two  capacitors  which 
exchange  heat  by  radiation 
GE0W)  -  surface  area 

F0FX  -  variable  name  given  to  the  table  listing  the  surface 
emissivity 

X  -  variable  name  given  to  the  table  listing  the  corresponding 
temperature 

10.  CALLWRITE  (N,  TIME1,  TIME2, . .  .TDiEN) 

N  -  number  of  times  at  which  printout  is  desired 

TIME1’ . TIMEN  *  specific  times  at  which  printout  is  desired 

11.  CALLASTEP 

This  subroutine  requires  no  calling  sequence. 

12.  CALLABL  (INDEX1 ,  INDEX2 , . . . . INDEXN) 

INDEX1,  INDEX2 , . . . . INDEXN  -  variable  names  assigned  in  control 

program  to  each  initial  surface  element 
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SECTION  4 
SAftruK  PROBLEM 

From  the  foregoing  discussions  of  the  program  structure  and  subroutine 
arguments,  the  ’•eader  should  now  have  a  general  idea  of  how  to  construct  a 
F0RTRAN  control  program  to  obtain  solutions  to  subliming  ablation  problems. 
However,  in  order  to  overcome  the  difficulties  inherent  in  the  word  descrip¬ 
tion  of  any  computer  program,  a  samnle  problem  will  now  be  solved  to  demonstrate 
precisely  the  use  of  the  computer  program. 

Pictured  in  Figure  4  is  a  sketch  of  the  area  to  be  analyzed.  The  area 
encompasses  the  juncture  between  a  thin  skin  and  a  thick  ring  both  of  which 
are  protected  by  .25  inch  thick  Teflon  ablator  material.  It  is  desired  to  find 
the  transient  temperatures  and  ablation  rates  in  the  ablator  and  the  transient 
temperatures  in  the  steel  understructure  for  a  Mach  5  low  altitude  flight.  The 
F0RTRAN  coding  needed  to  construct  the  main  control  program  will  be  described  in 
detail  and,  in  fact,  the  coding  sheets  themselves  are  included  in  this  section. 

No  attempt  will  be  made  to  discuss  or  justify  the  values  chosen  for  various 
parameters  (transition  Reynolds  number,  ablation  temperature,  etc.)  since  the  sole 
purpose  of  this  sample  problem  is  to  provide  instruction  in  use  of  the  program. 

It  is  usual  practice  to  divide  the  area  to  be  analyzed  into  a  network  of 
thermal  capacitors  before  the  coding  of  the  F0RTPAN  sheets  begins.  Previous 
experience  in  geometry  make-up  is  quite  helpful  in  order  to  make  a  good  selection 
of  element  size.  Naturally,  for  the  utmost  accuracy  the  smaller  the  element  size 
the  better.  However,  since  the  running  time  on  the  computer  is  increased  with 
decreasing  element  size,  a  compromise  between  accurac''  and  cost  must  many  times 
be  made.  It  is  sometimes  true  that  the  increased  accur-uy  due  to  a  smaller 
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break-up  of  the  geometry  is  not  measurable.  It  is  in  cases  such  as  this  that 
experience  is  helpful  in  determining  the  best  element  size.  The  geometry  break¬ 
up  chosen  for  the  sample  problem  is  shown  in  Figure  5.  With  the  brtsK-up  as 
shown,  the  solution  to  this  sample  problem  required  2.5  minutes  of  IBM  7094 
computer  time.  Had  the  break-up  been  twice  as  fine  as  that  shown,  the  running 
time  on  the  computer  would  not  have  increased  appreciably  since  most  of  the 
computer  time  was  used  in  converting  the  F0KTRAN  control  program  to  machine 
language  and  not  in  computation.  In  subdividing  the  solid,  it  is  generally 
advisable  to  hand-calculate  a  typical  stability  time  step  to  be  sure  that  the 
smallest  remains  greater  than  0.05  second.  Experience  shows  that  time  steps 
less  than  0.05  second  can  cause  excessive  consumption  of  computer  time.  Once 
the  element  sizes  have  been  chosen,  a  number  must  be  assigned  to  each  element 
and  the  volumes  of  each  element  and  the  area-to- length  ratios  between  adjacent 
elements  must  be  computed.  These  are  then  put  aside  to  be  included  later  in  the 
argument  of  the  CALIACAP,  CALLAO0MCN,  and  CALLO0N  statements.  A  detailed  dis¬ 
cussion  of  the  coding  sheets  now  follows.  The  list  of  units  which  must  be 
sed  with  the  i  lputs  is  given  in  Table  I  (Page  52). 

The  J0B  card  appearing  on  Page  1  of  the  coding  sheets  is  an  I.D.  card  whose 
form  is  peculiar  to  APL.  Other  systems  may  require  other  I.D.  cards  and  this 
one  is  shown  only  for  illustrative  purpose;.  The  XEQ  card,  however,  is  a  control 
card  and  must  appear  at  the  head  of  any  F0RTRAN  program.  The  SAMPIE  ABLATION 
PROBLEM  card  is  a  comment  card  which  is  not  processed  by  the  machine  and  is  used 
only  to  clarify  the  coding.  Conment  cards  will  be  used  liberally  throughout  the 
program.  The  first  of  the  two  DIMENSI0N  statements  describes  the  variables  listed 
in  each  particular  program.  Although  this  appears  at  the  head  of  the  program,  it 
is  usually  not  completed  until  the  rest  of  the  program  has  been  written  and  the 
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names  and  number  of  spaces  needed  for  each  variable  ere  known.  The  second 
DIMENSI0N  statement  must  appear  in  every  Job  in  exactly  the  form  as  shown  on 
the  first  coding  sheet.  Likewise,  the  C0MM0N  statement  must  appear  in  every 
job  exactly  as  shown  on  the  coding  sheet.  Since  the  second  DIMENSI0N  and  the 
O0MM0N  statements  are  never  to  be  changed,  they  may  be  placed  on  file  and  need 
not  be  re-coded  for  each  job. 

The  actual  coding  of  the  tabular  input  begins  on  the  second  coding  sheet. 

Here  the  Mach  number  versus  time  is  given  under  the  assigned  variable  names  FI 
and  F1V.  There  are  two  things  to  rote  in  this  table  which  are  true  for  each 
and  every  set  of  dependent-independent  variable  tables  written  in  the  control 
program.  First,  the  first  number  in  the  independent  variable  list  (F1V)  must  be 
the  number  of  items  in  the  dependent  variable  (FI)  table.  This  number  is  needed 
by  the  routine  which  makes  the  quadratic  interpolations  in  the  table.  The  second 
item  to  be  noted  is  how  the  points  describing  the  Mach  number  history  are  specified. 
In  this  particular  case,  a  linear  variation  of  Mach  number  for  5  seconds  up  to  a 
Mach  number  of  5  is  to  be  programmed.  After  5  seconds,  a  constant  Mach  number  of 
5  is  desired.  Anytime  such  a  sharp  break  occurs,  it  is  imperative  that  a 
number  of  points  extremely  .close  together  be  specified  in  the  area  of  this  break. 
This  is  necessary  because  quadratic  interpolation  is  automatically  used  between 
points  given  in  the  table.  If  enough  points  are  not  given  in  regions  where  there 
is  a  sudden  change  in  the  slope  of  the  curves,  completely  erroneous  interpolated 
results  can  be  obtained.  On  the  third  coding  sheet,  the  altitude  table  (E2)  is 
given.  No  new  corresponding  time  table  is  needed  since  the  altitude  values  were 
made  to  coincide  with  the  time  table  (FlV)  given  with  the  Mach  number.  Thus,  later 
the  Fl.V  table  will  be  used  with  both  FI  and  F2  in  the  argument  of  the  CALLTRAJ 
statement. 


! 
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Following  the  trajectory  data  the  ambient  atmospheric  conditions  are  listed. 
The  temperatuc and  pressure  are  specified  as  a  function  of  altitude.  !i»te  once 
again  that  the  coding  effort  is  reduced  by  making  one  altitude  table  serve  both 
the  temperature  and  pressure  tables. 

Following  the  ambient  pressure  table,  the  local  flow  conditions  (M/Mq  and 
P/Pq)  are  specified  as  a  function  of  free  stream  Mach  number.  It  might  be  noted 
here  that  unless  otherwise  stated,  the  variable  names  assigned  to  any  of  the 
tables  which  have  been  discussed  are  quite  arbitrary  as  long  as  floating  point 
names  (any  beginning  with  letters  other  than  I,  J,  K,  L,  M,  and  N)  are  used. 

For  example,  the  Mach  number  ratio  table  F5  could  equally  well  be  called  S5,  F27, 
or  T7. 

Following  the  local  pressure  ratio  table,  the  effective  space  temperature 
versus  altitude  is  listed  in  tabular  form  in  an  identical  fashion  to  the  other 
table?  discussed.  On  Page  6  of  the  coding  sheets,  the  emissivity  of  external 
surfaces  is  given.  In  this  case  it  was  desired  chat  a  constant  value  of  .8  be 
used  rather  tnan  a  temperature-dependent  set  of  values.  Note  that  since  this 
quantity  is  a  single  non-subscripted  constant  it  ne.,d  not  be  listed  in  the 
DIMENS I0N  statement. 

The  table  listing  the  flow  numbers  for  the  aerodynamic  heating  equations 
follows  next  and  is  the  only  table  which  is  the  same  length  for  every  problem. 
There  are  invariably  six  entries  in  the  tabic  and  each  of  these  entries  has  a 
particular  meaning.  The  first  entry  must  list  the  transition  Reynolds  number. 

The  second  and  third  entries  give  the  exponent  on  the  P.eynold3  number  in  the 
Nusselt  relation.  It  should  be  pointed  out  that,  within  the  routine,  the 
exponent  on  the  Prandtl  number  is  fixed  at  1/3.  The  foo'th  and  fifth  flow 
numbers  give  the  coefficients  of  the  Nusselt  number  relation  for  laminar  and 
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turbulent  flow  respectively.  The  values  assigned  to  the  second  through  fifth 
flow  numbers  in  the  sample  problem  apply  to  the  case  of  the  flat  plate.  The 
sixth  flow  number  is  set  equal  to  1  if  stagnation  point  heat  transfer  is  to  be 
computed.  At  all  other  times,  the  sixth  flow  number  must  be  set  equal  to  zero. 

The  next  table  lists  the  initial  volumes  of  all  the  potentially  ablating 
capacitors.  The  variable  name  used  in  this  list  is  not  arbitrary  but  must 
always  be  V0.  In  similar  fashion,  the  sequence  of  ablating  capacitors  in  the 
following  table  must  always  have  IASEQ  as  its  variable  name.  The  first  item  in 
this  list  must  specify  the  number  of  "columns"  of  ablating  capacitors.  For  the 
sample  problem  this  number  is  2  (see  Figure  5).  The  second  item  in  the  list  must 
specify  a  number  of  capacitors  in  the  first  column  (i.e.,  11).  Next,  the 
capacitor  number  assigned  to  each  of  the  elements  in  the  ablating  sequence  are 
listed  in  order.  Note  that  the  first  capacitor  in  the  understructure  (i.e.,  21) 

(see  Figure  5)  is  the  last  number  to  be  listed.  No  particular  number  sequence 
need  be  assigned  to  the  capacitors  when  formulating  the  geometry  break-up  as  any 
combination  of  numbers,  each  less  than  the  number  1000,  can  be  handled  by  the 
ablation  seque. ce.  The  fourteenth  entry  in  the  sample  problem  ablation  sequence 
is  the  number  of  capacitors  (i.e.,  11)  in  the  second  ablating  column.  The  fifteenth 
entry  starts  a  list  of  the  capacitors  in  the  second  column.  If  additional 
columns  were  present  the  cycle  would  be  repeated.  However,  the  total  length 
of  the  table  must  not  exceed  500  entries. 

Following  tin  ablatior  sequence  are  tables  giving  tb  thermal  properties  of 
the  two  materials,  Teflon  and  stainless  steel.  The  form  of  the  tables  is  identical 
to  those  listing  the  trajectory  and  atmospheric  properties-  Note  that  the  density, 
specific  heat,  and  thermal  conductivity  must  be  Jistc-d  individually.  In  the  case 
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of  stainless  steel  a  table  was  not  used  for  the  density  since  only  a  single 
value  was  kr.rwn.  This  value  is  given  on  Page  9  of  the  coding  sheets.  Constants 
may  of  course,  be  used  in  place  of  tables  at  any  time.  Following  the  listing 
of  the  thermal  properties,  the  ablative  properties  of  the  materials  are  given 
on  Page  11.  O.ice  again,  the  tables  are  in  the  usual  form  and  a  single  enthalpy 
difference  table  is  given  to  serve  both  the  HEFF  and  ablation  temperature,  tables. 

At  the  bottom  of  Page  11,  a  constant  equal  to  1  is  defined  under  the 
variable  name  of  F18,  This  constant  will  be  found  useful  later  in  the  argument 
of  some  of  the  subroutine  calling  statements.  A  final  item  specified  before 
the  CALL  statements  are  the  variable  names  assigned  to  the  initial  surface 
elements.  For  the  sample  problem,  element  1  has  been  assigned  the  fixed-point 
variable  name  INDX1  and  element  11  has  been  assigned  the  variable  name  INDX2. 

This  can  be  seen  on  Page  12. 

As  a  final  note  on  the  tables  it  should  be  mentioned  that  the  order  in  which 
all  the  above  tables  are  specified  is  not  ciitical.  They  may  be  placed  anywhere 
t etween  the  C0MM0N  statement  and  the  first  CALL  statement. 

The  remainder  of  the  program  consists  of  CALL  statements.  The  reader  will 
find  it  helpful  to  refer  periodically  to  Section  3  where  the  calling  sequence 
is  defined  for  each  of  the  routines.  Particular  attention  should  be  given  to  the 
way  variable  names  are  placed  in  the  calling  sequences  and  to  seeing  that  they 
do,  indeed,  satisfy  the  requirements  described  in  Section  3. 

The  first  statement  in  the  series  is  the  CALLASET  statement  (Page  12  of  the 
coding  sheets).  For  the  sample  problem,  CALLASET  has  been  instructed  to  make 
computations  from  0  to  11  seconds  and  set  all  capacitors  at  an  Initial  tempers- 
ture  of  519°R  with  the  exception  of  one  capacitor,  fu-aPered  1000,  whose 
temperature  is  set  at  600°R.  The  CALLTRAJ  statement  has  in  its  calling  sequence 
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the  variable  names  which  were  given  previously  to  the  Mach  number,  altitude, 
and  time  tablet  .  Note  that  this  statement  has  the  number  1  in  column  5.  It 
is  important  to  include  this  statement  number  since  the  program  is  instructed 
to  loop  back  to  statement  1  after  each  time  step.  The  CALLATM  statement  has  in 
its  calling  sequence  the  variable  names  given  to  the  ambient  temperature,  pres¬ 
sure,  and  altitude  tables.  The  CALLF0RALT  statement  lists  the  capacitor  number 
1000  which  was  chosen  as  the  capacitor  to  represenL  space  for  this  problem.  It 

also  lists  the  variable  names  given  to  the  temperature  of  space  and  to  the 

altitude  tables. 

The  CALLAAER0  statements  supply  information  to  compute  the  aerodynamic 
heating.  The  reader  should  refer  back  to  the  detailed  description  of  the  calling 
sequence  given  in  Section  3  and  to  the  accompanying  definitions  of  the  variable 

names.  In  the  arguments  of  the  CALLAAER0  statements  it  is  seen  that  surface 

areas  in  the  sample  problem  have  been  assumed  to  be  0.0833  square  feet  and  that 
the  reference  lengths  of  INDEX1  and  INDEX2  have  been  assumed  to  be  0.033  feet  and 
0,117  feet,  respectively. 

The  CALLAAP  statements  follow  next  on  Page  12  of  the  coding  forms.  Note 
that  there  is  one  statement  for  erch  capacitor  shown  in  Figure  5.  In  the 
calling  sequence  of  ACAP  the  hand-calculated  volume  of  each  capacitor  is  entered 
along  with  the  variable  names  assigned  to  the  density,  specific  heat,  and 
temperature  tables. 

On  coding  form  Page  13  are  listed  several  CALLC0N  statements  describing  the 
conduction  between  nonablative  capacitors.  CALLC0N  statements  are  used  in  place  of 
CALLAO0MCN  because  of  their  shorter  calling  sequence - 

Pie  CALLAC0MCN  statements  on  Page  14  identify  the  '’apacitor  pair  being 
joined  and  give  the  area-to-length  ratio  between  the  first  node  and  the  interface, 
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the  interface  area,  and  the  area-to-length  ratio  between  the  interface  and 

second  noce.  For  the  sample  problem,  it  was  assumed  that  the  interface  resis¬ 
ts 

tance  would  be  negligible.  Therefore,  a  ve-y  large  number,  .1  x  10' ,  was  used 

to  represent  the  interface  area.  The  interface  unit  conductance  was  then  assumed 

to  be  equal  to  1  and  was  inserted  into  the  calling  sequence  by  use  of  the 

previously  defined  variable  F18.  The  resulting  interface  conductance  is  so 
9 

large  (.1  x  10  )  that  the  interface  resistance  is  effectively  zero.  Note  that 
since  F18  was  a  constant  and  not  a  table,  there  was  no  independent  variable. 

A  zero  was  placed  in  the  calling  sequence  following  F18.  This  was  repeated  in 
subsequent  AO0MCN  sequences.  In  general,  if  a  constant  previously  assigned  a 
variable  name  is  to  be  entered  into  any  calling  sequence,  its  variable  name  is 
always  listed  followed  by  a  zero. 

On  Page  15  of  the  coding  forms  CALLRAD  statements  are  given;  one  for  INDXl 
and  the  other  for  INDX2.  In  this  problem,  1000  was  used  as  the  capacitor 
number  to  represen;  space  and  therefore  appears  in  the  calling  sequence  c  eac.. 

n 

statement.  A  surface  area  of  .0833  Ft“  is  given  followed  by  the  variable  name 
(F8)  assigned  ;o  the  constant  emissivity.  It  should  be  noted  once  again  that 
since  a  variable  named  constant  ic  used,  the  variable  name  is  followed  by  a  zero. 
The  CALLWRITE  statement  appears  on  Page  15  of  the  coding  sheets.  It  gives 
the  times  at  which  printouts  are  desired.  Care  must  be  taken  to  insure  that  the 
first  number  in  the  calling  sequence  is  an  accurate  count  of  the  number  of  times 
at  which  printout  is  desired.  One  should  recall  from  previous  discussions  thac 
printout  will  occur  automatically  at  crucial  times  other  than  those  lasted  in 
the  argument  of  the  CALLWRITE  statement.  In  order  to  allow  for  possible  longer 
running  times  future  printouts  were  called  for  up  t:  30  seconds  even  though  the 
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sample  problem  was  set  to  run  only  11  seconds,  lue  CALLASTEP  statement  following 
on  Page  16  of  the  coding  sheets  requires  no  calling  sequence  and,  therefore,  no 
discussion.  The  f^nal  call  statement  is  the  CALLABL  statement.  This  calling 
sequence  must  contain  the  variable  names  given  to  the  initial  surface  elements. 
For  this  sample  problem  the  names  INDX1  and  INDX2  were  assigned  to  the  initial 
surface  elements.  The  next  statement,  the  G0  T0  1  statement,  will  cause  the  pro¬ 
gram  to  go  back  to  the  CALLTRAJ  statement  after  each  time  step  and  proceed  with 
new  computations  for  the  new  flight  time.  The  program  is  completed  with  an  END 
card  which  is  required  in  all  F0RTRAN  jobs. 

After  the  data  on  the  coding  sheets  1  through  16  have  been  converted  to 
information  on  punched  cards,  the  resulting  control  program  deck  is  combined 
with  the  deck  of  binary  subroutines  for  running  on  the  computer.  The  decks  of 
binary  subroutines  are  available  upon  request  from  the  authors. 

Figure  6  presents  a  section  of  printout  taken  from  the  sample  problem  at  a 
flight  time  of  3.6  seconds.  For  the  most  part,  the  printout  is  self-explanatory 
but  a  few  words  of  clarification  will  be  given. 

The  units  of  the  variables  in  the  printout  are: 

Time  -  Seconds 

Altitude  -  Feet 


Pressure  -  Pounds  per  square  foot 

Temperature  -  °F 

H  prime  -  Btu-(Ft)a/Ft  • Sec*0F 

_  3 

Current  Volume  (VC)  -  Ft 

Effective  heat  of  ablation  (EFFHV)  -  Btu/lb 

Tne  thermal  capacitor  numbers  and  their  associated  temperatures  are  listed  if  the 
temperatures  have  changed  from  their  inicial  values.  Additional  contact  temperature. 
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are  given  on  both  sides  of  the  interface  between  those  capacitor  previously 
joined  by  CALLAC0MCN  statements.  The  local  flow  properties  are  listed  according 
to  an  ID  number  which  corresponds  to  the  ID  number  in  each  CALLAA.KR0  statement. 

Thus,  to  find  the  capacitor  numbers  to  which  each  set  of  flow  conditions  applies, 
one  must  refer  back  to  a  CALLAaER0  statement  containing  the  proper  ID  number. 

Figures  7  and  8  graphically  present  results  from  the  sample  problem.  Tn 
Figure  7  the  ablator  thickness  versus  time  for  the  two  ablating  "co±umn3"  is 
given.  The  ablation  rates  are  considerably  different  because  the  problem  was 
deliberately  set  up  with  boundary  layer  transition  occurring  between  the  two 
columns.  Thus,  v1ation  under  both  laminar  and  turbulent  flow  is  present.  Figure 
8  presents  temperature-time  results  for  the  external  surface  (at  the  location  of 
element  1)  and  the  steel  understructure.  It  is  seen  that  the  surface  temperature 
continues  to  rise  after  the  onset  of  ablation  at  3.4  seconds.  This  is  caused  by 
the  changing  enthalpy  difference  across  the  boundary  layer  and  the  fact  that  the 
ablation  temperature  is  a  function  of  the  enthalpy  difference. 

The  example  problem  chosen  is  an  extreme  case  which  demonstrates  the  capability 
of  the  computer  program  to  deal  with  a  substantial  differential  in  local  ablation 
rates  and  large  local  differences  in  the  thermal  mass  of  the  structure.  Therefore, 
it  should  be  re-emphasized  that  the  effect  of  body  contour  changes  on  vehicle 
aerodynamics  is  not  accounted  in  the  computer  program.  On  the  other  hand,  in  cases 
where  body  contour  changes  do  not  substantially  affect  vehicle  aerodynamics  the 
present  computer  program  can  be  used  effectively  to  integrate  the  thermal  analysis 
of  ablatively  protected  structures  having  complex  interconnected  conducting  paths. 

As  a  final  note  on  the  program,  Table  II  has  been  prepared  in  order  to  show 
the  upper  limit  on  program  capacity. 
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NOMENCLATURE 


a 


A 


C 


CP(s) 


a w 


W,i 


(g)w 


H 

H(s)o 

H(sg) 

H(sg)w 

H 

w 

H  . 
w,i 

* 

H 


Elemental  planar  surface,  located  in  conductive  median,  a.id  orientated 
in  direction  normal  to  direction  of  heat  flow  but  aligned  with  x- 
direction  on  ablation  element 

Elemental  planar  surface,  located  in  conductive  medium  and  orientated 
in  direction  normal  to  direction  of  heat  flow  but  aligned  with  z- 
direction  (in  or  out  of  exterior  surface)  of  ablation  element 

A  constant 

Capacitance  of  element  i 

Specific  heat  of  a  gas  at  constant  pressure 
Specific  heat  of  a  solid  material 

A  radiation  configuration  factor  based  on  geometry 

A  radiation  configuration  factor  based  on  emissivity 

Convective  heat  transfer  coefficient  based  on  enthalpy  difference 

Convective  heat  transfer  coefficient  based  on  enthalpy  difference  in 
the  absence  of  mass  release  from  surface 

Enthalpy  of  air  in  contact  with  an  adiabatic  wall  surface 

Enthalpy  of  air  in  contact  with  the  adiabatic  surface  of  element  i 

Enthalpy  of  gaseous  wall  material  at  wall  surface  temperature  and 
pressure  conditions 

Static  enthalpy  of  air  at  local-flow  conditions 
Enthalpy  of  solid  wall  material  at  initial  condition 

Enthalpy  of  sublimation 

Enthalpy  of  sublimation  for  the  wall  surface  temperature  condition 
Enthalpy  of  air  in  contact  with  a  non-adiabatic  wall  surface 
Enthalpy  of  air  in  contact  with  the  non-adiab- ric  surface,  qf  elpmpnt  -f 

Reference  enthalpy  of  air  defined  by  Equations  (3)  and  (4) 


1 
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ki  '  conductivity  of  element  i  at  its  temperature  condition 

k  ■  Thermal  conductivity  of  air 

K  -  Thermal  conductance 

"  Thermal  conductance  between  node  i  and  node  1 

Ki,il  *  T,)enDal  conductance  between  node  i  and  interface  between  element  i  and 

element  1 


i 


i,i3 


L 


i.il 


M 

o 

M 

m 

j 

n 

P 

P 

o 

Pr 


Elemental  length  of  conductive  path  directed  normal  to  adjacent 
aerodynamical ly  heated  wall  surface 


Elemental  length  of  conductive  path  directed  normal  to  adjacent 
aerodynamical!^  heated  wall  surface  and  running  between  node  i  and 
the  interface  between  element  i  and  element  3 


Elemental  length  of  conductive  path  directed  parallel 
aerodynamically  heated  wall  surface 


to  adjacent 


Elemental  length  of  conductive  path  directed  parallel  to  adjacent 
aerodynamically  heated  wall  surface  and  running  between  node  i  and  the 
interface  between  element  i  and  element  1 

Free-stream  Mach  number 


Mach  numbec  of  local  flow 
A  constant 


-  Mass  -low  rate  of  ablator  across  unit  area  of  the  control  surface 

-  A  constant 

-  Local  static  pressure 

-  Free-stream  static  pressure 

-  Prandtl  number 


-  Mon-ablative  convective  heat  flow  per  unit  time  and  unit  wall -surface 


-  Radiation  heat  flow  per  unit  time  and  unit  wall-surface  area 

-  Reynolds  number 

-  Absolute  temperature 
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T  -A  constant  used  in  Sutherland's  Viscosity  Law 

c 

Tq  -  Absolute  initial  temperature  of  solid  material 

T  -  Absolute  temperature  of  the  air  in  contact  w'th  an  adiabatic  wall 

AW 

“  surface 


T 

aw, 


i 


-  Absolute  temperature  of  the  air  in  contact  with  the  adiabatic  wall 
surface  of  element  i 


Tw  -  Absolute  temperature  of  a  non-adiabatic  wall  surface 

T  -  Absolute  temperature  of  effective  radiation  space 

space 


Aabl 


Absolute  temperature  of  an  ablating  surface 


^abl  i  "  Absolute  temperature  of  the  ablating  element  i 
T^  -  Present-time  absolute  temperature  of  element  i 

T*  -  Future-time  absolute  temperature  of  element  i 

k 

T  -  The  absolute  temperature  corresponding  to  reference  enthalpy  H  at 

1  atm. 

t,At  -  Time,  time  increment 

u  -  Local  velocity  external  to  boundary  layer 


V.  -  Presen. -time  volume  of  element  i 
1 

V+  -  Future-time  volume  of  element  i 

l 

x  -  Distance  from  leading  edge 


3jj  -  Transpiration  factor  based  on  enthalpy  difference 

k  * 

H  -  Air  absolute  viscosity  evaluated  at  T 

H  -  A  constant  in  Sutherland's  Viscosity  Law 

G 
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p  -  Air  density  evaluated  at  T 

a  -  Stefan-Boltzmann  Constant 
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TABLE  I 

UNITS  OF  COMPUTER  PROGRAM  INPUT  VARIABLES 


VARIABLE 
Altitude 
Temperature 
Pressure 
Volume 
Density 
Specific  Heat 
Thermal  Conductivity 

Modified  Effective  Heat  of  Ablatio:-  (HEFF) 
Area 


UNITS 

Ft 

°R 

Lbf/Ft2 

Ft3 

Lb  /Ft3 
m 

Btu/Lb  °R 
m 

Btu-Ft/Ft2 .hr • °R 
Btu/Lb 


Length 


Ft 
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APPEtPIX  A 

THE  MODIFIED  EFFECTIVE  HEAT  OF  ABLATION 

Figure  Al  has  been  prepared  in  order  to  show  two  typical  pyrolization 
processes,  and  a  typical  temperature  distribution  in  the  boundary  layer  and  ablator. 
Figure  Ala  illustrates  a  combined  sublimation-oxidation  process  while  Figure  Alb 
depicts  charring  >\«b  tot  bow;  .  With  reference  to  the  temperature  distribution.  Figure 
Ale,  it  can  be  assiaaed  either  that  the  ablator  has  been  suddenly  injected  into  a 
hot-gas  stream  and  is  experiencing  a  transient  surface  temperature  rise,  or  that 
ablation  has  steadied  out  such  that  the  temperature  profile  shown  is  of  a  fixed 
shape  and  magnitude  with  respect  to  the  moving  surface.  As  is  well  known,  one  of 
the  principal  features  of  ablation  is  the  reduction  la  inward  convective  heat  ft'.’.x 
through  the  boundary  layer  by  the  outward  flow  of  released  gas,  3his  release  of  gas 
into  the  boundary  layer  is  similar  in  its  effect  on  the  boundary  layer  to  the  release 
from  a  transpiratior  cooled  wall.  Therefore,  the  blocking  effectiveness  of  the 
release  is  commonly  characterized  by  a  transpiration  factor 

n 
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Figure  Al.  Tvo  Modes  of  Ablation  and  a  Typical  Temperature  Distribution 
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In  the  case  of  the  steady  ablation  process  the  fundamental  steady 
one -dimensional  energy  equation  on  t- e  control  volume  of  Figure  Ala,  which 
extends  deep  into  the  steadily  ablating  semi-infinite  body,  is 


=  H,  ,  -  H.  ,  +  p„(H  -  H  )  +  Mrad 

(g)w  (s)o  H  aw  w  — r~ 

m 


Effective 
Heat  of 
Ablation 


Enthalpy 
Rise  of 
Ablator 


Blocking  Net  Outward 
Term  Surface 

Radiation 


where  qQ  denotes  ordinary  unblocked  convective  heat  flux  to  a  wall  whose  surface 

temperature  is  equal  to  that  of  the  ablating  surface.  With  reference  to  equation 

(Al)  it  can  be  seen  that  a  purely  experimental  effective  heat  of  ablation  can  be 

obtained  in  the  laboratory  by  applying  a  given  level  of  (H  -  H  )  and  by  measuring 

aw  w 

qo  to  a  non-ablating  calorimeter  model,  whose  surface  temperature  is  equal  to 
the  ablator  surface  temperature,  and  by  measuring  mass  loss  rate  m  of  the  ablating 
model  in  the  same  stream.  The  experimentally  determined  quotient  qo/'m  will 
correctly  represent  the  right-hand  side  of  equation  (Al)  if 

(a)  Calorimeter  radiation  matches  that  of  the  ablator,  or  if  qQ  is 
corrected  in  the  case  of  a  mis -match. 

(b)  The  ablation  model  is  thick  enough  to  simulate  a  semi-infinite  body. 

(c)  The  run  is  steady  or  is  corrected  for  any  unsteady  periods. 

(d)  Ihe  supply  air  is  clean. 

(e)  Air  pressure  effects  on  degradation  chemistry  simulate  the  in-flight 


By  comparison  to  the  foregoing  direct  experimental  m.-hods,  the  indirect  experimental 


method  of  determining  qQ/m  i.wolves  selecting  a  level  ^ r  HaW-  Hw  and  suitable 


T 


TN»  Mr*  HnUw 
<miM  hwm  UHunn 


experimental  values  ol  other  quantities  on  the  right-hand  side  of  Equation  (Al). 
In  the  particular  case  of  a  subliming  ablator,  ’oe  calculated 

from  the  relation 


H(g)v  H(s)o  *  H(sg)*r  *  Jj,  Cp(«) 


c  /vvdT 


provided  the  ablation  wall  surface  temperature  is  known.  Hie  enthalpy-based 

transpiration  factor  S0  is  available  in  the  literature  (5)  covering  the i laminar 
n 

and  turbulent  boundary  layers.  However,  calculation  of  the  radiation  term  qQ/ni 
is  not  straight-forward  because  it  involves  the  unknown  ablation  rate  m.  The 
alternate  expression  for  qo/m  due  to  Adams  (5) 


H.  .  +  -  H  ) 

(s)o  H  aw  v' 


results  from  a  rearrangement  of  the  terms  of  Iquation  (Al).  Still,  however,  the 
calculation  of  .  ./q  is  not  straight-forward  because  it  involves  the  unknown 

heating  rate  q^.  No  attempt  will  be  made  to  resolve  this  problem  here.  For 
further  discussion,  see  Reference  6. 

In  the  present  ablation  program  an  equivalent  surface-ablation  energy 
property  was  needed.  It  was  defined  as  the  sum  of  the  sublimation  enthalpy 
change  and  the  blocking  term.  By  substituting  Equation  (A2)  into  Equation  (Al) 
one  finds  this  sum  to  be  related  to  the  effective  heat  of  ablation  in  the  following 


I  %  r>T w 

H(sg>w  +  ^H^aw  "  *(  m  jj 

'  0 


semi -infinite 
test  body 
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Since  the  righ.-hand  side  of  equation  <A4)  is  a  modified  form  of  the  ordinary 
effective  heat  of  ablation  qo/n  it  is  called  the  modified  effeccive  heat  of 
ablation.  When  obtaining  the  modified  effective  heat  of  ablation  from  the 
literature  one  must  exercise  good  judgement  with  regard  to  factors  (a)  -  (e) 


cited  above. 
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